function [A,M] = FEMMatrix(V,T,kappa,dt)
%
% pre-calculate the finite element matrix using linear finite element
% on quadrilatteral
%
ne = size(T,1);
idx = reshape([1 2 3 4]'*ones(1,4), 16, 1);
jdx = reshape(ones(4,1)*[1 2 3 4], 16, 1);

is = zeros(ne*16,1);
js = is;
As = is;
Ms = is;
l = 1;
for i = 1:ne
    nd = T(i,:);
    vtx = V(nd,:);
    
    is(l:l+15) = nd(idx);
    js(l:l+15) = nd(jdx);
    
    [mass,stiff] = elementMat(vtx,kappa,dt);
    
    Ms(l:l+15) = reshape(mass,16,1);
    As(l:l+15) = reshape(mass,16,1) + dt*kappa*reshape(stiff,16,1);
    
    l = l + 16;
end

ndof = max(max(T));
M = sparse(is, js, Ms, ndof, ndof);
A = sparse(is, js, As, ndof, ndof);
end


% element matrix of linear finite element on quadrilaterral with vertex vtx
% use linear fem on triangle element -- 1 to 2 for quadratural
function [Mass,Stiff] = elementMat(vtx,kappa,dt)
% for stiff
R_11 = [2 -2 -1  1; -2  2  1 -1; -1  1 2 -2;  1 -1 -2  2]/6;
R_12 = [1  1 -1 -1; -1 -1  1  1; -1 -1 1  1;  1  1 -1 -1]/4;
R_22 = [2  1 -1 -2;  1  2 -2 -1; -1 -2 2  1; -2 -1  1  2]/6;
% for mass 
R_00 = [4 2 1 2; 2 4 2 1; 1 2 4 2; 2 1 2 4]/36;

F = inv([vtx(2,:)-vtx(1,:); vtx(4,:)-vtx(1,:)]);

E = F'*F;

Mass = R_00/det(F);

Stiff = R_00 + kappa*dt*(E(1,1)*R_11 + E(1,2)*R_12  + E(2,1)*R_12' + E(2,2)*R_22)*det(F); 


end

